A graphical theory of competition on spatial resource gradients 
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Resource competition is a fundamental interaction in natural communities. However little is 
known about competition in spatial environments where organisms are able to regulate resource 
distributions. Here, we analyze the competition of two consumers for two resources in a one- 
dimensional habitat in which the resources are supplied from opposite sides. We show that the 
success of an invading species crucially depends on the slope of the resource gradients shaped 
by the resident. Our analysis reveals that parameter combinations which lead to coexistence in 
a uniform environment may favor alternative stable states in a spatial system, and vice versa. 
Furthermore, differences in growth rate, mortality or dispersal abilities allow a consumer to coexist 
stationarily with - or even outcompete - a competitor with lower resource requirements. Applying 
our theory to a phytoplankton model, we explain shifts in the community structure that are induced 
by environmental changes. 
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INTRODUCTION 

Competition for limiting resources is one of the most 
important species interactions in ecology and has long 
been considered as a major driver for shaping the spa- 
tial structure of communities and limiting species rich- 
ness. A graphical theory of resource competition was ad- 
vanced by MacArthur (1972), Leon & Tumpson (1975) 
and Tilman (1980, 1982) and revolves around the as- 
sumption that consumers that reduce limiting resources 
to the lowest level will exclude all other competitors (the 
i?*-rule). Being confirmed experimentally (Miller et al. 
2005, but see also Wilson et al. 2007), this theory pro- 
vides a fundamental framework for interpreting the re- 
lationship between organisms and their shared resources 
in a uniform environment (Grover 1997; Chase & Lei- 
bold 2003). Extensions of competition theory to spatial 
settings (Chesson 2000a, 2000b; Klausmeier & Tilman 
2002; Amarasekare 2003) can be grouped into different 
classes. In spatially homogeneous environments, coexis- 
tence can be mediated by trade-offs in life history pa- 
rameters (Levins & Culver 1971; Levin 1974; Kneitel & 
Chase 2004). In heterogeneous environments regional co- 
existence can be mediated by spatial segregation (Tilman 
1982; Gross and Cardinale 2007), whereas local coex- 
istence may be promoted via source-sink effects (Mou- 
quet & Loreau 2003; Leibold et al. 2004) or positive 
correlations between dispersal and competitive abilities 
(Abrams & Wilson 2004). However, previous theory has 
not thoroughly explored resource competition in contin- 
uous spatially variable habitats. 

In a system with resource flows, spatial heterogene- 



ity can be created by biotic interactions, as local re- 
source consumption can modify resource levels over larger 
ranges (Huston & DeAngelis 1994), possibly establish- 
ing resource gradients over the full extent of the habitat. 
The analysis of competition in such systems is mathe- 
matically challenging (Grover 1997), as not only the re- 
source availability, but also the size of favorable areas 
(Skellam 1951), dispersal rates (Abrams & Wilson 2004) 
and many other factors become crucial for the survival 
of a population (Ryabov & Blasius 2008). Despite exten- 
sive theoretical studies on unstirred chemostats (Hsu & 
Waltman 1993; Smith & Wahman 1995; Wu et al. 2004), 
persistence and competition in streams (Speirs & Gurney 
2001; Lutscher et al. 2007), vegetation patterns in water- 
limited habitats (Klausmeier 1999; von Hardenberg et al. 
2001), well-mixed systems with a light gradient (Huisman 
& Weissing 1994; Diehl 2002), and non- uniform phyto- 
plankton systems (Dutkiewicz et al. 2009; Yoshiyama et 
al. 2009; Ryabov et al. 2010), translating competition 
theory to extended systems, where species are able to 
shape resource distributions, still remains a great chal- 
lenge. 

In this study we develop a general framework for ana- 
lyzing the competition between two consumers for two 
limiting resources in a spatially continuous habitat in 
which resource distributions are regulated by biotic in- 
teractions. In extension of the i?*-rule for uniform sys- 
tems, we introduce the notion of an invasion threshold 
in a spatial system as the maximal resource requirement 
for a consumer to invade in the presence of a resident 
species. In this way, the outcome of competition can be 
interpreted graphically in the resource plane, by compar- 
ing the location of the invader's critical resources with 
respect to the invasion threshold lines. We derive ana- 
lytic expressions for the invasibility conditions and show 
that they can be related to how resources are spatially 
reduced as a result of resource-use by the resident. 

Using this approach we find that environmental hetero- 
geneity increases the likelihood of coexistence. Thereby 
we identify two forms of stationary coexistence in spa- 
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tially variable habitats. First, coexistence can be me- 
diated by a trade-ofF in resource requirements. This is 
characterized by spatial segregation of the two species 
and we find that parameter combinations which allow 
coexistence in a uniform system favor alternative sta- 
ble states in a spatially extended system and vice versa. 
Second, coexistence can arise from positive correlations 
of growth rates and resource requirements (gleaner- 
opportunist trade-ofF), characterized by a lack of spatial 
segregation. 

To illustrate these ideas, we investigate two phyto- 
plankton species competing for light and a nutrient in 
a water column. We identify two distinct regimes of 
bistability in the competition outcome and show that our 
approach provides a powerful basis for projecting the in- 
fluence of environmental changes on the community com- 
position. 



MODELLING FRAMEWORK AND ANALYSIS 
Resource competition in uniform systems 

As a background for the following discussion, we briefly 
review the theory of resource competition in a uniform 
environment (MacArthur 1972; Leon and Tumpson 1975; 
Tilman 1980, 1982; Grover 1997; Chase & Leibold 2003). 
Consider two consumer species i = 1 , 2 growing on two 
essential resources N and /. The dynamics of the pop- 
ulation density Pi of each species is determined by the 
difference between growth /i^ (N, I) and mortality rrii 



P, = {^M{NJ)~m,)P,, 
where the growth rate in general can be written as 



lJLi{NJ) = ^„ 
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Here Umax^i is the maximal growth rate of species i, /i 
is a monotonically increasing function of both arguments 
with upper bound liuiiqj^oo = 1, and the half- 

saturation constants i?Ar,i and Hj i define the species' 
resource adaptation. In particular, this functional form 
includes von Liebig's law of minimum, where the limita- 
tion of growth follows the Monod kinetics 
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N + H, 
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(3) 



In a uniform environment, levels of system resources 
can be represented as a point in a two-dimensional re- 
source plane (Fig. la). In this plane, states of balance 
between growth and mortality, fJ-i{N, I) = mj, determine 
the zero net growth isoclines (ZNGIs), which for each 
species divide the resource plane into areas of positive 
and negative population net growth. Assuming eqn[3]for 
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Figure 1: Invasion analysis in a uniform system (top) and 
a spatially extended system with opposing resource gradi- 
ents (bottom), (a) Equilibrium configuration of the resident 
(species 1, green). The system state point at equilibrium E 
is determined by the intersection of the consumption vector 
CV, starting at the resource supply point S, with the zero 
net growth isocline ZNGI (dashed line, see text), (b) Critical 
resource values of a successful invader (species 2, red) should 
be located below the invasion threshold T (blue line). The 
intersection of the invader's ZNGI with T is shown as cross- 
hatching, (c) In a spatially extended system, combinations of 
resources at different spatial coordinates give rise to a system 
state curve SSC (solid line, favorable range is marked in bold), 
(d) Invasion threshold T (blue solid line) and first order ap- 
proximation (blue dashed line) with slope 71 = ci,i/c]v,i in 
the log-log plot (eanll2[). (c) and (d) show results of numer- 
ical simulations in the phytoplankton model, eqns 3-[6] (see 
Appendix S5 for model parameters). 



the growth rate, we can calculate the critical values, /* = 
Hi,i rrii/ {fimax,i - mi) and N* ^ H^,, m,/ {^imax,^ - rrii), 
and the ZNGI takes the form of two orthogonal lines 
(dashed lines in Fig. 1). 

The resource concentrations in the absence of con- 
sumers define the coordinates of a resource supply point 
S. A growing population will deplete resources, shift- 
ing the system state in the direction of the consumption 
vector CV , until the system state point hits the popula- 
tion's ZNGI at the equilibrium point E — {Ni,Ii). This 
point defines the resource levels at equilibrium, which are 
shaped by the population, and sets the conditions for the 
invasion of other species. 

Assume that resident species 1 has attained equilib- 
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rium. Then invasion of a second species 2 is possible 
if the equilibrium resource concentrations shaped by the 
resident are sufficient to allow positive growth of the in- 
vader. Thus, in resource space a part of the invader's 
ZNGI must be located below the equilibrium point E 
(Fig.lTJ)). For the specific form of eqn 3, this condition is 
fulfilled for all critical resource values of the invader that 
are located in the rectangular area defined by < Ii 
and N2 < Ni . This line defines the invasion threshold T 
for species 2 (blue line in Fig. lb). 

Combining the invasion analysis for each of the two 
species yields the well-known outcomes of resource com- 
petition in a uniform environment (Tilman 1982). With- 
out a trade-off in the use of the two resources (i.e. the two 
ZNGIs do not intersect), the species with the lowest re- 
source requirements wins. In the presence of a trade-off, 
competition can lead to stable coexistence, competitive 
exclusion or bistability (i.e. alternative stable states with 
outcomes depending on initial conditions), depending on 
whether or not each species has a greater impact on the 
resource that most limits its own growth (see Fig. SI). 



Resource distributions in a spatial system 

The classical approach implicitly assumes the presence 
of boundaries which confine the system, between which 
organisms and resources are uniformly distributed. In 
the following we aim to translate these ideas to an open 
system where the favorable area of a species is confined 
by resource availability, and where organisms and re- 
sources are not distributed uniformly. Consider a one- 
dimensional environment, where the population densities 
Pi{z,t) depend on a spatial coordinate z. A spatial ana- 
log to eqn 1 can be written in the form 



dt 



[ix,{N, I)-m,)P, + D, 



(4) 



where Di characterizes the diffusive dispersal ability of 
population i. Again the growth rate fJ,i{N, I) follows the 
general form of eqn [2l but the resource concentrations 
N{z,t) and I{z,t) are functions of z, so that the two 
populations are indirectly coupled through the spatial 
profile of their shared resources. 

As a case study, we examine a two-species 
phytoplankton-nutrient model in an incompletely mixed 
water column (Radach & Maier-Reimer 1975, Jamart et 
al. 1977, Klausmeier & Litchman 2001, Huisman et al. 
2006, Ryabov et al. 2010). In the model, eqn 4 describes 
the density Pi{z, t) of phytoplankton species i at depth z 
and time t, and is complemented by an equation for the 
nutrient N{z, t) 



D 



92 AT 

dz^ 



(5) 



and an equation for the light intensity I{z,t), which de- 
scribes the absorption of light by water and phytoplank- 
ton (Kirk 1994) 



I{z) = hn exp 



-Kbgz- / y2hP,{^,t)d^ 



(6) 



In these equations, parameter D obtains the role of the 
turbulent diffusivity, is the nutrient content of a phy- 
toplankton cell, lin is the incident light intensity, Kbg is 
the water turbidity and ki is the attenuation coefficient 
of phytoplankton cells. As boundary conditions we as- 
sumed impenetrable borders at the surface and at the 
bottom for the phytoplankton biomass and an impene- 
trable surface and a constant concentration Nb at the 
bottom for the nutrient (for model parameters see Table 
SI). 

The model describes a situation in which the two re- 
sources, light and nutrient, are supplied from opposite 
sides of a spatial habitat. This is typical in the ocean 
where light is supplied from above and many macronu- 
trients from below (Fig. 2a). As shown in Fig. 2b, this 
gives rise to characteristic resource distributions with in- 
verse gradients, where population growth is maximal at 
an intermediate position at which both resources are suf- 
ficiently available (see Fig. S6 for further typical model 
solutions). 

Within the favorable range of the population, the light 
and nutrient distributions in Figs. 2a and 2b can be well 
approximated as straight lines on a logarithmic scale. 
This means that resource distributions at equilibrium de- 
cay (or grow) exponentially in space 



N,iz) 



(7) 



This exponential dependence can be derived analytically 
in the limit of low mortality (see Appendix SI), and it 
is typical for many other ecosystems (see e.g. Fig. S5). 
Eqn 7 implies that a population i in monoculture is 
"shading" resources by a constant percentage 



CN,: 



1 dN,{z) 
Ni dz ■ 



c/,. 



1 dhjz) 
h dz 



(8) 



Here, the logarithmic resource gradients CM,i and c/_j 
measure the infiuence of the population on the resource 
distribution. 

Translated into the resource plane, the resource dis- 
tributions of the spatially extended system can be com- 
pactly represented as a parametric curve, (iV(z), /(z)), 
(solid line in Fig. Ic). This 'system state curve' (SSC) 
naturally extends the notion of a system state point E 
of a homogeneous system to a spatially continuous envi- 
ronment (Ryabov et al. 2010). As shown in Fig. Ic, due 
to the spatial coupling and source-sink effects the system 
state curve at equilibrium does not settle at the ZNGI. 
Instead, it intersects the ZNGI at two points which mark 
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Figure 2: Typical distributions of light (gray), nutrient (black) and chlorophyll/phytoplankton (green) concentrations in a 
water column from (a) field measurements (HOT program, Station ALOHA, cruise 114; see also Karl 2010) and (b) numerical 
simulation of eqns [SHS] Dashed lines in (a) and (b) indicate exponential fit of the resource distributions, (c) Net growth rates, 
fii{N, I) — m, of a resident (species 1, green) and an invader (species 2, red), derived from the simulated distribution of resources 
(solid lines) and from their exponential approximation (eqn 7, dashed lines). Horizontal dashed lines in (b) and (c) mark the 
favorable range of species 1. Vertical shifts Azn and Azi, related to nutrient and light limitation (eqn 10), are also shown. 
Model parameters as in Fig. 1. 



the boundaries of positive net growth. This inner part of 
the system state curve (marked in bold in Fig. Ic) cor- 
responds to the resident's favorable range (horizontally 
dashed lines in Fig. 2b, c). Given the exponential depen- 
dence (eqn 7) , within the favorable range the system state 
curve approaches a straight line with slope 7i = ci^i/cN,i 
in the resource plane with logarithmic axes (see Fig. Ic 
in the model or Fig. S4c for field data). 



Invasion analysis in a spatial system 

In analogy to competition theory in a uniform system, 
we may now ask whether it is possible to predict the out- 
come of two-species competition from knowledge about 
equilibrium resource distribution which has been shaped 
by each species alone in the absence of the other species. 
Assuming that we know the resident's system state curve 
at equilibrium (iVi(z), /i(z)), what can we say about the 
success of an invader (red color in Fig. 2c)? 

To address this problem we define the invasion thresh- 
old T in the resource plane as the set of critical resources 
of all invaders of small initial density, which have zero to- 
tal growth in the presence of the resident (blue solid line 
in Fig. Id, obtained by invasion analysis in numerical sim- 
ulations with 7000 parameter combinations). In this for- 
mulation, the success of an invader depends on whether 
or not its ZNGI intersects with the invasion threshold 



(cross-hatching in Fig. 1). Thus, invasion analysis can 
be performed graphically by comparing the location of 
the invader's critical resources with the invasion thresh- 
old. 

By comparison of Figs, lb and Id, the invasion thresh- 
old has a different shape in a uniform or a spatially 
continuous system. It is difficult to predict the inva- 
sion threshold from general principles. Exact conditions 
can be expressed in terms of eigenvalues of a correspond- 
ing boundary value problem (Hsu & Waltman 1993; Ap- 
pendix S3) , which unfortunately cannot be solved in gen- 
eral. In the following we show that, with the assumption 
of exponential resource distributions (eqn 7), it is pos- 
sible to calculate the invasion threshold for species with 
similar resource requirements. 



Calculation of the invasion threshold 

Consider first the simplest case where resident (species 
1) and invader (species 2) differ only in their half satura- 
tion constants, Hjsi^i and Hj i, but are otherwise identical 
(Mmaa;,i = t^^maxa, "ii = m2, £>i = D2) ■ As the luvadcr 
density is assumed to be small, and its influence on the 
resource distributions can thus be neglected, the possibil- 
ity of invasion depends entirely on the invader's growth 
rate in the equilibrium resource distribution shaped by 
the resident, pL2{Ni{z)^Ii{z)). The mathematical iden- 
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tity 



J_gC2 ^ J_gC(z-A2 



1 

with Az = - In — , (9) 
c iJi 



H2 Hi 

shows that division of an exponential function by differ- 
ent (half-saturation) constants Hi and H2 is equivalent 
to a shift Az in position along the z-axis. Therefore, 
assuming that the growth rate follows eqn [5] and the re- 
source distributions are exponential (eqn [7]), it is possible 
to express the growth rate, fJ-2{z), of the invader through 
that of the resident 



(10) 



H2 (Ni (z), /i (z)) = ^1 (iVi (z - Azjv) , h (z - Az/)) 

Here the shifts Az^r — c^^j^ ln_ffjv,2/^Ar,i and Az/ = 
—cj \ \nHi^2/Hi^i have opposite signs, since the resource 
gradients are inverse. 

Introducing^ a newposition z' = z — Azn, we find that 
)i2iz) = fii{Niiz'),h{z' + A)), with A = AzAT - Az/. 
The change of variables z ^> z' corresponds to a shift 
along the z-axis (see Fig. 2c) and has no effect on the 
conditions for survival as long as boundary effects can 
be neglected. Thus, the difference in growth between in- 
vader and resident depends only on the value and sign 
of the single parameter A. If A = 0, the distinct re- 
source requirements of the species just result in a paral- 
lel translation of the growth rate profile. Since the res- 
ident species at equilibrium has zero total growth, the 
same holds for the invader, and the population of the in- 
vader cannot establish. By contrast, if A > 0, we obtain 
I{z' -I- A) < /(z') since /(z) decays with z. This leads 
to negative net growth of the invader since the function 
IJ,{N, I) increases monotonically with its arguments. Us- 
ing the same arguments, A < gives rise to positive 
growth of the invader. 

Thus, we can formulate the invasibility criterion in the 
form A < 0. The parameter A describes the invasion 
threshold and can be interpreted in terms of resource 
gradients and relative competitive abilities 



A 



CNA 



In- 



H^ 



N.2 



Hn- 



1 

CIA 



In 



Hij. 
Hi- 



inverse log of relative 
log gradient competitive 
in N abilities for N 



inverse log of relative 
log gradient competitive 
in / abilities for I 



(11) 

If for instance, both species have equal resource require- 
ments for iV, His[,2 = Hpf^i, invasion is only possible, 
A < 0, if the invading species is a better competitor for 
/, i.e. ln(ff/,2/ff/,i) < 0. 

If the limitation of growth follows von Liebig's law 
(eqn [3]) , the spatial profiles for light and nutrient limi- 
tation are shifted independently and the net shift A — 
Azn — Azi measures the difference in the size of the res- 
ident and invader's favorable habitats. This is depicted 
in Fig. 2c for a situation with Azat < Azi, where the 
narrowing of the favorable range due to higher N re- 
quirements {N2 > N*) is less than the widening due to 



better / adaptation (/^ < ). Therefore the net growth 
rate of the invader is greater than that of the resident, al- 
lowing for successful invasion (see Fig. 2c). Note that in 
the example shown in Figs. 1 and 2, the invader, species 
2, is a good light competitor, /| < /j, but needs higher 
nutrient concentrations, N2 > N^. Therefore it will have 
maximal production at a deeper depth (Fig. 2c). 

Since the maximal growth and mortality rates are as- 
sumed to be equal for both species, the ratio of half- 
saturation constants equals the ratio of critical resource 
values (e.g.. Hi 2/ Hi i = I2/I1), and the invasibility cri- 
terion reads 



A = 



CN,1 N* 



1 T* 
— ln^<0. 

c/,1 It 



(12) 



This expression has a straightforward geometrical inter- 
pretation: a species can invade the spatial habitat if its 
critical resources {N2, /|) are located below a straight 
line, In / — In /j" = -7i(lniV — InA^^*), with a slope of 
absolute value 71 = ci,i/cn,i passing through the point 
{Nl, II) in the double- logarithmic resource plane (blue 
dashed line in Fig. [TJl). Thus, we have derived a first- 
order approximation for the invasion threshold T. 



Competition outcome 

In a similar way it is possible to analyze the invasion 
potential of species 1 in a system with resident species 
2. Assume without loss of generality that iVj > and 
I2 < Ii, and denote by jcr the slope (taken with op- 
posite sign) of a straight line passing through the two 
critical resource points, 7cr = -(ln/|//j')/(ln AflZ-^i)- 
Then, combining eqn [T^] and its counterpart for species 
2, we obtain four different outcomes of spatial resource 
competition (Fig. 3a). If the critical resource values for 
each species are located below the invasion threshold of 
its competitor, 71 < jcr and I/72 < l/7cr, both species 
can invade the monoculture of the other, leading to co- 
existence. In the opposite case, if the critical resource 
values for every species lie above the invasion threshold 
of its competitor, 71 > jcr and I/72 > l/7cr7 neither 
of the two species can invade, leading to bistability. Fi- 
nally, one species may be a superior competitor if only 
this species can grow in the presence of its competitor. 

This analytic conclusion is confirmed by numerical sim- 
ulations of the phytoplankton model eqns 3-6 in a large 
parameter range (Fig. 3b). The figure demonstrates that 
in spite of the high complexity and nonlinearity of the 
model, the relation between the slope of the invasion 
threshold lines, 74, and the location of the critical re- 
source values gives a good prediction of the competition 
outcome. 
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Figure 3: The four different cases of spatial resource competition depending on the ratio 7i of logarithmic resource gradients: 
(71 < 7cr, 1/72 < l/7cr) coexistence; (71 > -jcr, I/72 > l/7cr) alternative stable states; (71 > 7cr, I/72 < l/7cr) species 
1 wins; (71 < jcr, I/72 > l/7cr) species 2 wins, (a) Schematic representation. Insets sketch invasion threshold lines for 
resident species 1 (green) and invader species 2 (red), and critical resource values (dots) in the resource plane, (b) Numerical 
verification in the phytoplankton model, based on more than 1000 simulations, with parameter values randomly chosen from 
a wide range (critical resources (A',*, I*) and 7cr remain constant, see Appendix S5). The color scale shows the logarithm of 
relative species abundance \nBi/B2 at equilibrium (total biomass Bi = J^^ Pi{z)dz). Upward-pointing triangles indicate the 
simulation outcome when species 1 is the resident and species 2 the invader, downward-pointing triangles show the opposite 
case; their color is different only if the competition outcome is bistable. 



(a) 



Successful 
invasion 



(b) 




N 



Figure 4: (a) Interspecific differences in jj-max, rn, or D yield 
a shift 5i2 of the invasion threshold in the normal direction. 
Compared with Fig. Id, the invader (species 2, red) can es- 
tablish, even though it has higher resource requirements than 
the resident (species 1, green), (b) Coexistence due to a 
gleaner-opportunist trade-off, where both invasion thresholds 
are shifted in opposite directions. Species 1 (gleaner) with 
lower resource requirements can coexist at equilibrium with 
species 2 (opportunist), which has higher resource require- 
ments but higher fj,max, or lower m or D. 



Differences in growth, mortality and dispersal 

In general, the locations of the invasion thresholds may 
depend on trait differences between the resident and in- 



vader, such as maximal growth rate, mortality and dis- 
persal rates. As shown in Appendix S3, the generalized 
form of the invasibility criterion reads 



CAT,! N* 



— Inf <Ai2 (13) 



where A12 is a complicated function of the trait values 
and resource gradients, but does not depend on the crit- 
ical resource levels. If all trait differences vanish, we ob- 
tain A12 = and recover the previous result, eqn 12. The 
geometrical representation of the new term A12 is a shift 
of the invasion threshold line toward larger (A12 > 0) or 
smaller (A12 < 0) resource requirements in the normal 

direction by 612 = Ai2Cn,ici^i/ \J c\[ i + i- Thus, the 
slope of the invasion threshold lines remain unchanged, 
but their location is shifted. Fig. 4a shows an exam- 
ple with a positive shift of the invasion threshold, where 
species 2 can invade the system despite the fact that it 
has higher requirements for both resources than the res- 
ident species 1. 

For an illustration of how specific trait values can 
change the sign of ^12, consider the particular case that 
species 2 has higher values of both growth rate, ^imax,2 = 
P^imax,!, and mortality, 7712 = /3mi, with j3 > \. Since 
the critical resource values N2 > Ni and /| > /j* are 
independent of /3, this scaling would not affect the com- 
petition outcome in a uniform environment and species 
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(a) 



1st wins 



Bistability: 




2nd wins 
1st or 2nd 



Coexistence 
1st or coex. 



Figure 5: Competition analysis in the phytoplankton model. 
Maximal growth rates and mortalities are (a) identical, /3 — 
1, and (b) different, 13 — 2. The competition outcome (in 
color coding) of species 1 (critical resource values indicated as 
circle) is shown in dependence of the critical resource values of 
species 2. Solid lines show invasion thresholds obtained from 
numerical simulation of eqns |3Hn] dashed lines in (a) show 
their analytical approximations (Appendix S2). Parameter 
combinations {N2 , 12) in the area below the solid green line 
allow invasion of species 2 in the presence of species 1. In 
contrast, species 1 can invade in the presence of species 2 with 
critical resource values taken from the area above the solid 
red line. The intersections of these lines define the ranges of 
coexistence, alternative states, and competitive exclusion. 



2 would not be able to invade the system. However, in a 
heterogeneous system the invasion threshold (green line 
in Fig. 4b) is shifted by S12 > (see Appendix S3). Simi- 
lar analysis, using rescaling by a factor of shows that 
the threshold for invasion of species 1 in a monoculture 
of species 2 is shifted in the opposite direction (red line in 
Fig. 4b). This mechanism can lead to stable coexistence 
of a species having high growth rate and resource require- 
ments with another species having low growth rate and 
resource requirements (gleaner-opportunist trade-off). 



Application to a phytoplankton model 

Applying this approach to a phytoplankton commu- 
nity, the slopes 7^ of the invasion lines can be obtained 
by numerical simulations (Fig. S2), analytic estimations 
(Appendix SI), or in field enclosure experiments (Jager 
et al. 2008). Tracing the dependence of 7^ as a function 
of system parameters allows us to project shifts in the 
species composition. Consider, for instance, the effects of 
increasing the nutrient resource concentration at the bot- 
tom of the water column, Nb ■ More nutrients at the bot- 
tom give rise to a larger biomass of the resident species, 
which in turn yields a steeper nutrient decay within the 
production layer, i.e. a reduction of ji (see Fig. S2c). 



In the resource plane this effectively leads to a counter- 
clockwise rotation of the invasion threshold line, and 
graphical analysis reveals that by increasing Nb the com- 
munity composition shifts from dominance of the best nu- 
trient competitor, through a regime of coexistence, to the 
prevalence of the best light competitor (Fig. S3). This 
result might explain observed negative correlations of the 
abundance of high light-adapted species with nutrient 
concentrations (Johnson et al. 2006). 

Similarly we can study the influence of consumption 
rates. In a well-mixed water column the slope of con- 
sumption vectors is given by Ti — [ki zb) / cti (Huisman 
& Weissing, 1994; Diehl 2002). Thus, the relative re- 
source consumption scales with the ratio of the light at- 
tenuation coefficient to the cell nutrient content, ki/ai. 
In our model, according to Fig. [31 two species coexist 
if 71 < 7cr < 72- Since the ratio of resource gradients 
7i also grows with ki/ai (see Fig. S2a and eqn S6), we 
observe coexistence if the best nutrient competitor (here 
species 1) has a small fci/ai ratio and species 2 a large 
k2/a2 ratio. The reverse situation (large fci/ai and small 
k2/a2) leads to bistability (Fig. [5^). Thus, we obtain a 
rule which is diametrically opposed to that for uniform 
systems: coexistence arises if each species reduces its 
least limiting resource more strongly in relation to the 
other, whereas stronger reduction of the most limiting 
resource favors bistability (compare Figs. [3] and SI). 

Fig. [5^ shows the competition outcome in dependence 
of the critical resources of species 2 for the case that 
ki/ai > fc2/a2, be. species 1 shades relatively more light 
than species 2. In accordance with our rule, if species 2, 
which has a smaller light attenuation coefficient, is more 
strongly limited by light {1 2 > /j*), it may coexist with 
species 1 (upper blue area), whereas stronger nutrient 
limitation of species 2 {N2 > N^) can lead to bistabil- 
ity (white area). However, if the difference in nutrient 
limitation is sufficiently large, then the two species can 
again coexist (lower blue area. Appendix S2). The model 
further exhibits a regime of bistability between coexis- 
tence and single-species dominance (gray area), where 
the invasion of species 1 does not exclude species 2, but 
species 2 cannot invade if species 1 has already estab- 
lished (see also the vertical profiles in Fig. S6). If the 
two species differ in their maximal growth and mortality 
rates, (3 > 1 (Fig. [5]d), all bifurcation lines are shifted 
towards larger critical resource values. Then, if species 
2 is parameterized to have higher resource requirements 
(i.e., critial resource levels located above the black dashed 
lines in Fig. \Sjp) , it can have all possible competition out- 
comes, even without a trade-off in resource requirements. 
The corresponding spatial profiles in Fig. S6 reveal two 
fundamentally different types of coexistence in a spatial 
system, mediated either by a trade-off in resource re- 
quirements (characterized by spatial segregation of the 
two species) or by a gleaner-opportunist trade-off (lack 
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of spatial segregation). 

DISCUSSION 

The extension of resource competition theory to spa- 
tially variable environments remains a major challenge 
for community ecology. Usual approaches to include en- 
vironmental heterogeneity consider a range of resource 
supply points, assuming that the system evolves indepen- 
dently for each supply point (Tilman 1982), while other 
studies have concentrated on patch occupancy mod- 
els (Levins 1979) or mctacommunity models of discrete 
patches that are connected by dispersal (Levin 1974; 
Mouquet & Loreau 2003; Abrams & Wilson 2004; Gross 
& Cardinalc 2007). However, links between local patches 
might not only affect the dispersal of species, but also 
the spatial flow of matter and resources, leading to meta- 
ecosystcm dynamics (Loreau et al. 2003) which require a 
distinct theoretical approach (Huston & DeAngelis 1994, 
Smith & Waltman 1995, Wu et al. 2004). 

Our study takes a significant step in this direction by 
exploring competition in a spatially continuous environ- 
ment where species interact indirectly by modifying re- 
source availability beyond their local neighborhood. We 
study the critical resource requirements of a successful in- 
vader in such a system and find that the invasion thresh- 
old can be approximated by a straight line on a double- 
logarithmic scale, with a slope that is determined by the 
ratio of logarithmic resource gradients. Thus, the like- 
lihood of invasion in a two-species community can be 
understood from the analysis of equilibrium resource dis- 
tributions in a monoculture of the resident. 

Using this approach, we clarify the bases of coexis- 
tence in spatially variable habitats for a model that has 
widespread and important applications (see Table 1). 
Thereby, we synthesize two previous sets of results con- 
cerning coexistence: i) uniform habitat theory, showing 
that coexistence can arise from trade-offs in the use of 
two resources (Tilman 1980, 1982), and ii) spatially vari- 
able theory showing that coexistence can arise from a 
growth-based trade-off (Hsu & Waltman 1993; Smith & 
Waltman 1995; Wu et al. 2004). We show that in a spa- 
tially variable habitat with two resources, both forms of 
coexistence can occur. 

Assuming that two species trade off in their resource 
adaptation, in a uniform system both species can coexist 
if each species mostly consumes its most limiting resource 
(Leon & Tumpson 1975). In contrast, in a spatial envi- 
ronment the species coexist if they each, in monoculture, 
create a resource distribution with a relatively smaller 
gradient (c/ or c^r) of their most limiting resource. In 
other words, invasion is possible if the resident species 
does not "shade" its most limiting resource too much. As 
a consequence, parameter combinations which allow co- 
existence in a uniform system can lead to alternative sta- 



ble states in a spatially extended system and vice versa. 

These striking differences between resource competi- 
tion in uniform and spatially variable systems can be 
explained by the observation that in a uniform en- 
vironment, resources are equally available everywhere, 
whereas in a spatially extended system consumers com- 
pete for locally available resources. Thus in a uniform 
system, a species can suppress invasions by reducing 
those resources for which it is the best competitor. In a 
spatial setting the same strategy may not work because 
a local reduction of vital resources does not exclude in- 
vasions elsewhere in the habitat. Consider for instance 
the situation shown in Fig. 2c. Here the resident, species 
1, is the best nutrient competitor, so that in a uniform 
environment invasions are prevented if species 1 reduces 
nutrients more strongly than light (i.e. a small value of 
ki/ai, or Fi, in Fig. SI). In a spatial situation, how- 
ever, this strategy (small value of fci/ai) will promote 
invasion of species 2 (see Fig. 3), because species 1 can 
only shade nutrients vertically above its favorable range, 
whereas species 2 invades at a greater depth. There- 
fore, to prevent invasions, species 1 shoiild reduce light 
as much as possible (large value of fci/ai) to deteriorate 
growth conditions in the deep layers. 

Applying our theory to a spatial phytoplankton model, 
we identify two distinct regions of bistability, of either 
alternative states for each species or bistability between 
coexistence and a monoculture. Bistability has been the- 
oretically described in phytoplankton communities with 
differently mixed layers (Yoshiyama et al. 2009, Ryabov 
et al. 2010). Here we identify bistability in a system 
with uniform diffusivity, however, unlike in well-mixed 
systems, the range of bistability is smaller than that of 
coexistence (see Fig. 5). 

Finally, in a spatial system the slope of the invasion 
threshold changes with the biomass of the resident (see 
Appendix SI), because species with higher density have a 
stronger influence on resource distributions. In this way, 
competitive interactions are intricately linked to produc- 
tivity, which potentially can give rise to novel causal re- 
lations between productivity and diversity (Gross & Car- 
dinal, 2007; Hillebrand & Matthiessen, 2009). 

Extensions of the theory 

An essential part of our analysis is based on the as- 
sumption of exponential decay of resources within the 
favorable range of the resident (eqn 7) . This assumption 
makes spatial invasion analysis analytically tractable, as 
it allows us to express the growth rate of an invader by 
that of a resident species, thereby circumventing the com- 
plicated problem of actually having to calculate spatial 
density proflles. Exponential resource distributions cor- 
relate with field data and can be observed in aquatic 
systems (Figs. 2a and S4; Kirk 1994; Karl & Letelier 
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Uniform system 


Spatial system with resource gradient 


Invasion threshold 


simple form defined by resource values, 
on the ZNGI 


curved shape, defined by gradients and resource availability, 
location depends on differences in fimax, m, or D 


Coexistence with 
resource trade-off 


each species mostly consumes 
its most limiting resource 


each species mostly consumes its least limiting resource or 
essential difference in resource requirements 

— > spatial profiles with separation 


Coexistence without 
resource trade-off 


not possible 


species with higher resource requirements (opportunist) 
has larger ^max, or smaller m or D 

— s- spatial profiles without separation 


Bistability 


equally likely as coexistence, 

each species mostly consumes 
its least limiting resource 


less likely than coexistence, 

each species mostly consumes its most limiting resource, 
possible also without resource trade-off 



Table I: Comparison of competition outcome in uniform and spatially extended systems with inverse resource gradients 



2008), but also appear in other systems, such as marine 
sediments (Fig. S5). 

With increasing distance from the favorable range 
one should expect deviations from exponential profiles 
(Fig. 2b). Invasion thresholds then become curved 
(Fig. Id) and our analytic theory only provides a first- 
order approximation for competition between biologically 
closely related species (but note that our general ap- 
proach of investigating invasion thresholds does not have 
this limitation). In this sense our theory complements 
the approach developed by Yoshiyama et al. (2009) for 
competition between sufficiently different species. 

For simplicity we assumed that the favorable ranges 
of both species are far from the system borders, so that 
boundary effects on the species survival are negligible. A 
general approach should include both forms of the inva- 
sion threshold, presented in Fig.[T|D and[T}l, as limit cases, 
where the range of positive growth is confined either by 
the boundaries of a well mixed system or by the resource 
gradients. Our numerical simulations in the phytoplank- 
ton model (not shown) demonstrate that alternative sta- 
ble states replace regimes of coexistence or vice versa, as 
the favorable range moves from the interior of the water 
column to the surface. 

The analysis can be extended to take other system pa- 
rameters into account. For example, in Appendix S4 
we show how sinking or fioating can effectively be in- 
terpreted as a change of mortality rates, thus also influ- 
encing the position of invasion thresholds. To simplify 
settings we considered competition between two species 
under equilibrium conditions. The approach could be 
generalized to include temporal changes in habitat con- 
ditions. Resource fluctuations shift SSCs together with 
invasion thresholds and can provide different time niches 
for r and K strategists (Grover 1990, 1991). A further 
interesting perspective would be to consider additional 
species. The success of a third competitor can be ana- 
lyzed by including another invasion threshold in the re- 
source plane. However, analysis of three-species com- 
petition is more intricate because it involves the study 
of invasion thresholds for three possible pairs of resident 



species. In preliminary investigations we have found that 
the combination of gleaner-opportunist and resource lim- 
itation trade-offs may promote stationary coexistence of 
three species on two resources. 

The method developed in this manuscript is indepen- 
dent of the nature of the abiotic resources involved. Sit- 
uations with spatial profiles of two inverse chemical re- 
source gradients arise naturally whenever resources are 
provided from two opposite sides, such as interfaces be- 
tween different environments (D'Hondt et al. 2004). 
Moreover, our theory can easily be extended to include 
spatial gradients of other vital factors, independently of 
whether or not they are actually infiuenced or consumed 
by the population. In fact, such a situation is obtained 
as a special case of our phytoplankton model by setting 
the light attenuation coefficient ki = 0, so that light in- 
tensity decays exponentially with depth in the water col- 
umn, independently of the biomass distributions. Sim- 
ilar spatial variation of growth factors or stress agents 
is abundant along environmental gradients and in tran- 
sition zones between different habitat types. Thus, the 
general approach that we develop should apply to a wide 
spectrum of systems, ranging from flow reactors, stream 
ecosystems and transport-limited vegetation systems to 
marine sediments and many others. 
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SUPPORTING INFORMATION 



SI. Calculation of the logarithmic resource gradients 

In this section we provide analytic estimates for the hght and nutrient gradients resuhing from a single species 
population in the phytoplankton model, eqns 3-6. We calculate resource distributions in the production layer, defined 
by the critical depths and Zj where either the nutrient or the light intensity reaches a critical value, N{z'^) = N* 
' ' , ' and I{zj) = I* (see dashed lines in Figs. 2b and c). We further assume small mortality m ^ iJ-max, so that the 

O ' '^^^^^ ''^"^ ^"^"^^^^ phytoplankton density has only small 

variation within the favorable layer, i.e. P{z) « Pq for z"^ < z < z*j. 

• From eqn 6 within the favorable range the logarithmic gradient of the light intensity equals 
(D ■ 

c,^y-^^K,, + kPo, (SI) 

and the light distribution can be written 
PLh ■ I{z)^ Pe^'^-^'-'K (S2) 

d 
cr 



Regarding the nutrient profile we assume that for small mortality the critical nutrient concentrations are small, 
N* <C H]y, so that the growth rate close to the critical point N w N* can be linearized fiN{N) w iJ-max^ff^ ■ 
Substituting this expression into eqn 5 we obtain at equilibrium 

^ N dN 

>^ ' -o:^i,nax-r;—Po + D-r^ = ■ 

■ Hn dz'' 

^, 

A solution to this equation, that is monotonically increasing with depth, is given by 
g '. N{z) = Ar*e="(^-^«) ^ (S3) 

with the logarithmic gradient 



To complete our calculation the phytoplankton density in the production layer is estimated as Pq ~ ■§ where B is 
the total biomass and the width of the production layer w = |z^ — z|| only weakly depends on the model parameters 
(Beckmann and Hense 2007). If the biomass maximum is located far from the surface, the value of B can be estimated 
in two limiting cases (Klausmeier and Litchman 2001) 



iiPok:^Kbg 



B = { (S5) 

>{Nb~N" ' ' 



where Zm = Zb - ln{Irn/I*)/Kbg. 



These equations can be further simphfied because typically Nb ^ N* . Substituting eqns [S5] into eqns [ 
we obtain the ratio of resource gradients 'y = ci / cn (see also Fig. 

ln(/,;„//*) 



2 

I and El 



bg 



Hhn/I* 



if Pofc > Kb 



k DH 



N 



7= \ 



(S6) 



k DNb 



a wmZ 



M 



Nb 



if Pnk < Kb 



Hff wmZ 



M 



S2. Bifurcation lines in Fig. 5 

In this section we derive analytic estimates for combinations of critical resource values that allow mutual invasibility 
of a fixed "reference" species 1 with a variable "test" species 2 (bifurcation lines in Fig. 5a main text). Both species 
differ in the values of consumption rates, ki and a^, have identical maximal growth rate fJ-max and mortality m, but 
while the critical resource values N* and Ii of the reference species are fixed, N2 and I2 are taken from a wide range. 

The critical resource values of all test species that can invade in the presence of the reference species are determined 
by the invasibility criterion, eqn 12 

1 N* ^ T* 

.ln^ + ^h4<0. (S7) 



The border of this region is shown as a green dashed line in Fig. 5a. Similar, we obtain the critical resource values of 
the test species, which allow the invasion of the reference species when the test species has reached equilibrium 

1 /V* 1 A* 

+ — 1^4<0. (S8) 



CN,2 C/,2 /I 

In this case however, to plot the bifurcation line, we need to consider the dependence of cm,2 and C7.2 on the critical 
resource values, N2 and /|, of the test species (see. eqns ISH IS4I and [S5|) . For these aims we pick first one special 
"hybrid" test species, that has the same critical resource values Nl, as the reference species, but consumption 
rates, 02 and k2- We denote the logarithmic resources gradients of this species as c'^ 2 and c'j 2- 

To proceed we assume that the critical resource values are small: assuming N* <^ Nb we approximate Nb — N.2 ~ 
Nb, and assuming /* ^ /;„ we obtain 

T T / T* — T* 

m — — « m C ' 



r* T* \ T. 

Thus, in both limits, Kbg ^ Pq/c and Kbg <C Pok, in eqn [S5] the total biomass B, and from eqn ISII also the value of 
c/^2j is independent of the critical resource values iV|, /| of the test species. This means that the logarithmic light 
gradients are identical for all test species 



c/,2 = c'j 2 



3 



In a similar way we estimate the logarithmic nutrient gradients, but from eqn[S4]the value of cn,2 depends on the 



half-saturation constant of the test species, cn.2 ^ \/^/ Hn,2- Then the logarithmic nutrient gradient of an arbitrary 
test species can be expressed in terms of the logarithmic nutrient gradient of the "hybrid" species 



CN,2 — Cjv,2 

and using the identity ffjv,i/^^jv,2 — NI/N2 we obtain 



H 



N,2 



Cn,2 



-N,2 





Therefore, inequality IS8I takes the form 

7^1n§>0, (S9) 
'1 '-^1,2 ^1 

where the parameters 2 '^'i 2 calculated for the "hybrid" species having the critical resource values, and 
Ii, and the consumption rates, k2 and a2. Note that this expression holds in both limits when Kbg ^ Pok and when 
Kbg Pak, implying that it might provide a good approximation in general. The border of this region is shown as 
red dashed line in Fig. 5a. 



S3. Invasion thresholds and principal eigenvalues 



The dynamics of a spatially extended population P{z, t) is determined by a reaction-diffusion equation, eqn 4, with 
certain boundary conditions. The solution to this equation can be presented in the form 

00 

P(z,t) = ^e^"VnW , (SIO) 

n=l 

where ipn and A„ are the eigenfunctions and eigenvalues of the equation 

(/i(z)-m-A„)Vn + ^^^ = (Sll) 

with the same boundary conditions for the functions i}jn{z). The eigenvalue A* with the largest real part, the so-called 
principal eigenvalue, is of special interest for determining the dynamics of t) because the term with this eigenvalue 
dominates the sum ISlOl at large times (see e.g., Ryabov & Blasius 2008). The sign of the real part of the principal 
eigenvalue determines if the population will grow, and so is able to persist A* > 0, or if it declines exponentially on 
the whole habitat. A* < 0. Note, that A* depends on the spatial resource profiles and therefore may change in time 
with depletion of resources by a growing population. If a species has attained equilibrium, it has zero population net 
growth and its principal eigenvalue equals zero. A* = 0. 

The principal eigenvalue can be used for invasion analysis in a spatially extended system. We can define a family 
of invaders which for small initial density in the presence of the resident species have zero net growth, and thus a zero 
principal eigenvalue A*„^ = 0. In the space of the invader's critical resource values, this family defines the invasion 
threshold (in general, a hypersurface) , separating successful and unsuccessful invaders. In our setting, and given that 
resident and invader have identical maximal growth rate, /^max,2 = Mmax,i, and mortality, 7712 — mi, this corresponds 
to the condition A — Cjy\ \\iHi^^2/ Hn,i + Cj \ lni//_2/i?/,i = (see eqn 11). 
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Now consider an invader which has the same half saturation constants as the resident species, but different maximal 
growth rate, /imax.2, mortality, TO2, or dispersal ability, D2. It can be shown (Stakgold, 2000) that the eigenvalue 
spectrum, and in particular the principal eigenvalue, of this species shifts in the positive direction with an increase 
of the growth rate difference imax.2 — ^J■max.l, and in the negative direction with an increase in the difference of 
destructive factors, such as dispersal rates D2 — Di or mortalities m2 — mi. Thus, in general, the invader's principal 
eigenvalue will be positive or negative, even though the constant remains A = (since the half saturation constants 
remain unchanged). Assume without loss of generality that the invader has a positive total growth rate, A*„^ > 0. 
This means that it should be possible to invade the system with even higher resource requirements. Thus, there 
should be a 'border-line' invader with the same values of ^max,2, and D2, but larger half saturation constants so 
that its principal eigenvalue equals zero, A*„^ — 0. As the resource requirements of this second invader are higher, its 
A value is positive. This value, which we denote by AI2, defines the difference in the favorable ranges (according to 
eqn 11) between this species and the resident. Thus all species with maximal growth rate /^max.2, dispersal rate D2 
and mortality m2, and whose half-saturation constants satisfy the inequality 

In -— + — In -— < A12 (S12) 

can invade the system. As the boundary problem can not be solved in general, we cannot calculate the value of 
A*2 but we can make some predictions about its sign. Similar to A*„^, this value increases with {firnax,2 ~ l^max.i), 
decreases with m2 — mi and D2 — Di, equals zero when all these differences vanish, and changes sign together with 
A* 

inv ■ 

Assuming Monod limitation of growth by two essential resources (eqn 3) we can rewrite eqn S12 in terms of the 
critical resources values 

1 N* ^ T* 
-^ln^ + -Lln^<Ai2, (S13) 

CN,1 C7,i li 

where 

A12 = At2 + f — + — ) In -^""^^i"^^ . (S14) 

VCJV,1 Cl^iJ (Mmax,2 - ™2)™l 

Eqns. S13 and S14 define the invasion threshold for the invader, species 2, as a straight line with slope 71 = Cn,i/ci^i 
in double logarithmic resource space. The location of the invasion threshold depends in a complex way on the resource 
gradients and on the differences in maximal growth rate, /^max.i, mortality, rrii, or dispersal rates, Di. In the resource 
plane, the invasion threshold can be located above (A12 > 0) or below (A12 < 0) the critical resource values of the 
resident, (A^j*,/]*), see Fig 4a. 

Consider a special case when the growth and mortality rates of species 2 are rescaled, //max. 2 = /?Mmax,i and 
TO2 = /3mi. Then the second term in eqn IS14I vanishes, however A*2 can be nonzero because both Hmax.2 — l^max,i 
and 7712 — mi are positive. Substituting /imax,2 and m2 into eqn 4 and dividing it by /3, we obtain that in equilibrium 
the distribution of species 2 should satisfy the equation 

Did^P2 

fJ.i[z)P2 - miP2 + — — TTV = ■ 

Therefore these changes can be interpreted as a reduction of the diffusivity D2 — Di/ (3. Since A12 decreases with 
(Z?2 — Di) the shift of the invasibility threshold is positive {/S.\2 > 0) if /3 > 1 and negative (A*2 < 0) if /? < 1. 
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S4. Extension to systems with sinking or advection 

The dynamics of a population living in a unidirectional flow (e.g. sinking, floating, etc.) can be described by the 
reaction-diffusion-advection equation (Murray 2002, Ryabov & Blasius 2008) 

dP dP d^P 

_.,„M-,„,P-.._ + i,_, ,S15) 

where v is the flow or sinking velocity, which is positive when the flow is directed towards larger z values. Substituting 
P — Pexp {vz/2D) into eqn ISf 5[ we obtain 

dP ~ { v'^\ ~ d^P 

Note that dtP has the same sign as dtP, so that both functions grow and decline simultaneously. Furthermore, they 
have the same boundary conditions P{0) = P{Zb) = 0. Egn. ISf6l describes a population P{z,t) in a system without 
flow but with an effective mortality 

m' = m + . (Sf 7) 

4D 

Therefore, the presence of an advective flow which washes out a population from a favorable range can be interpreted 
as an effective increase of the mortality by the value /AD. 

This higher mortality entails an increase of resource requirements. Assuming eqn 3 for the growth rate, we obtain 
new critical resource values 

777' 77?' 

/r = i//,. — ^ , ivr = — '—^ . (SI8) 

Since m' > m, the new critical resources values /'* and N'* in the presence of a positive flow are larger than in a 
system without flow. This, in turn shifts the zero net growth isoclines and invasion thresholds towards higher resource 
values in the resource plane. 
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S5. Model details and parameters used for figures 

Numerical scheme The model was integrated using a backward difference method, based on the finite volume 
scheme (Pham Thi et al. 2005). For the numerical solution we have discretized all variables on a grid which consisted 
of 400 points. Diffusion terms were approximated by a second order central discretization scheme and integration 
was made via the trapezoidal rule. The resulting system of ordinary differential equations was solved by the CVODE 



package (http://www.netlib.org/odel. Further we verified that the results remain unchanged if we double the number 
of points in the grid. See Ryabov et al. (2010) for further details on numerical procedures. 

Boundary conditions As boundary conditions we assumed impenetrable borders at the surface and at the bottom 
for the phytoplankton biomass and an impenetrable surface and a constant concentration Nb at the bottom for the 
nutrient. 

Initial conditions We used a linear initial nutrient distribution changing from at the surface to Nb at the bottom. 
For both phytoplankton species we assigned a uniform distribution of small initial density. We also investigated the 
influence of different initial conditions, however we did not find any new solutions beside the ones described. 
To perform the invasion analysis, in the numerics the growth of the invader species was suppressed during the first 
10000 simulation days, to make sure that an equilibrium distribution of the resident species was established. Then 
we simulated the system for a duration of further 40000 days to obtain the final competition outcome. To test 
for bistability, this simulation was repeated by a second simulation in which the roles of invader and resident were 
exchanged. 

Parameter values Default species and model parameters are listed in Tabic [ST] 

Figure Id. Calculation of the invasion threshold is based on 7000 simulations of invasion by species with different 
half saturation constants Hjs! and Hi. However residents and invaders have the same /Zmax, w, and D, so that the 
invasion threshold is not shifted with respect to the resident {N*,I*) values. 

Figure 2. Half-saturation constants of the invader, i?}"^ = 8 jimol photons m^^ and i/]^^ =0.09 mmol nutrient m^''. 
Figure 3b. The figure shows the results of more than 1000 simulations where parameter values of environmental 
conditions K^g, D, Nb) and values of consumption rates, ki and a^, were randomly chosen from a uniform 

distribution in a wide range (see Table SI). Specific parameters, which do not vary, correspond to species 1 and 2 
in Table [SH In particular, the half-saturation constants, maximal growth and mortality rates were fixed, so that the 
critical values, {N* , I*), and jcr remain constant. To determine 71 and 72 numerically we performed simulations 
for each species in monoculture. Typically a change in environmental parameters leads to a simultaneous increase or 
decrease of 71 and 72 . As a result all data points group along one diagonal. 

In the case of competitive dominance (where one species wins and the other is excluded) the ratio p = InBi/ B2 
practically approaches positive or negative infinity after our standardized simulation period of 40000 days (see above). 
To visualize intermediate values of p close to the coexistence region, the color scale was truncated at absolute values 
of I/O I > 7, i.e. log-ratios above 7 where set to p = 7 and log-ratios below -7 where set to p = —7. 
To judge bistability, simulations were repeated two times, with either species 1 taken as resident and species 2 as 
invader (upward-pointing triangles) or species 2 as resident and species 1 as invader (downward-pointing triangles). 
Figure 5. Invasion thresholds are plotted on the basis of invasion analysis for 7000 combinations of species 1 and 
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test species 2. Parameter values taken as in Table ISll However, for species 2 half saturation constants were varied 
to achieve the chosen values of critical resources, and maximal growth rate and mortality taken as (a) ^max,2 = 0.04 
h-i and ma = 0.01 h^^ and (b) ^niax,2 = 0.08 h^^ and ma = 0.02 hr^ . 
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Figure SI: Two-species resource competition in a uniform system (compare to Fig. 3). The four typical cases for the 
competition outcome are shown, depending on the mutual slopes of consumption vectors Ti (values of critical resources and 
supply points remain fixed). For each case the zero- net-growth isoclines (ZNGI's, dashed lines) of species 1 (green) and 
species 2 (red) are plotted in the (A'^, /)-resource plane, together with the supply point (black circle) and the consumption 
vectors (green and red arrows). Each species in monoculture reduces resources to an equilibrium point, corresponding to the 
intersection of its consumption vector with its ZNGI. Invasion of the other species is possible if its critical resource values are 
located below this point (see Fig. lb). Define the critical slope Tcr by the slope of a straight line from the supply point to 
the intersection of the two ZNGIs (indicated as black dashed line). Species 1 can invade if the consumption vector of species 
2 has a slope less than the critical slope (l/r2 > l/Fcr, top panel), i.e. species 2 has relatively stronger influence on its most 
limiting resource A'^ than on resource I. In contrast, species 2 can invade if the slope of the consumption vector of species 1 is 
larger than the critical slope (Fi > Fcr, right-hand panel), i.e. species 1 has a relatively stronger influence on its most limiting 
resource /. 
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a D Ng 

Figure S2: Ratio of the logarithmic resource gradients 71 = c/,i /cm,i in the phytoplankton model as a function of several species 
specific and environmental parameters. Colorcode indicates the value of arctan 71 calculated numerically for a monoculture of 
species 1 in the phytoplankton model (see Table [STI for parameters which do not vary). Since 7^ equals the slope of the invasion 
line in the resource plane, this allows to project shifts in the species composition and to gain insight of how a change of one 
parameter can compensate for the influence of other parameters. For example, an increase of 7 corresponds to a clockwise 
rotation of the invasion threshold lines, thereby favoring the best nutrient competitor. 

(a) Dependence of 7 on species traits. Slope of the invasion threshold 7 grows with the ratio k/a. This dependence is also 
evident from both limits in egn IS6I Together with Fig. 3, this result conflrms our suggestion that two-species coexistence is 
more probable if each species relatively less consumes its most limiting resource (see main text). 

(b) Dependence of 7 on environmental conditions. Slope of the invasion threshold 7 grows with turbulent diffusivity D and 
background turbidity Kbg. This indicates that a turbid but weakly-mixed environment should result in the same species 
composition as stronger mixed but clean waters (while the total biomass will, of course, in general be different). 

(c) Dependence of 7 on resource supply. Increase of incident light intensity or of bottom nutrient concentration Nb decreases 
7 (thus leading to a counter-clockwise rotation of the invasion thresholds) and creates more favorable conditions for the best 
light competitor (compare to Fig. IS3|l. 
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Figure S3: Geometrical method to project the outcome of spatial resource competition in dependence of the ambient nutrient 
concentration, Nb- Top panel; shift in the competition outcome between species 1 (green) and species 2 (red) in the phyto- 
plankton model caused by an increase of the bottom nutrient concentration, Nb- Bottom panel illustrates this shift as the 
result of a counter-clockwise rotation of the invasion thresholds. The slopes, 71 and 72, of the invasion thresholds in the bottom 
panel were calculated numerically for a monoculture of species 1 and 2. To obtain a relatively small difference between 71 and 
72 we used ai = 8 x 10^^" mmol nutrient cell~^ and — 6x10"^" m^ cell~^. 
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Figure S4: Exponential resource distributions in a vertical water column. Shown are the profiles of photosynthetically active 
radiation (gray), nitrate (blue) and chlorophyll (green) measured during the HOT program (ht tp:/ /hahana.soest.hawaii.edii] 
station ALOHA, cruise 114; see also Karl 2010). Resource distributions plotted (a) on linear and (b) on logarithmic scales. As 
it is commonly supposed, the light intensity can be approximated by an exponential distribution (black dotted lines and 
values in (a) and (b)). In contrast, the nutrient profile exhibits two distinct regimes: In the area below the production layer 
the nutrient concentration can be well fit by a linear dependence {R^ — 0.91, black dashed line in (a)); however, within the 
production layer an exponential distribution gives a better description (i?^ — 0.99, black dashed line in (b)). We argue that 
the regime of exponential decay within the production layer is more crucial for competition, because in this area growth can 
be nutrient limited, whereas in the linear part of the nutrient profile growth is nutrient saturated. Interestingly, the interval of 
exponential decay ranges over three orders of nitrate concentrations, making this area crucial for the competition among a wide 
range of phytoplankton species which can be limited by this nutrient (Litchman et al. 2007). Note, that similar exponential 
distributions in this interval of depths have been observed in mean nitrate concentrations between 1989-2006 (Karl & Letelier 
2008). (c) System state curve (SSC), showing the distributions in the resource plane with logarithmic axes. The solid blue 
line shows actual cruise data for the light intensity and the nutrient. Since light intensities are not reported below a depth of 
140 m, the system state curve was extrapolated, assuming that the light intensity follows the same exponential decay up to 
300 m (blue dotted line). Dashed line shows a straight line fit (with R^ value) within the production layer. 
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Figure S5: Exponential distribution of resource concentrations in deep marine sediments. Inverse concentration profiles 
(note the logarithmic scale) result from biologically catalyzed reactions, which consume and release metabolites in a complex 
biogeochemical reaction network deep in the sediment column. 

(a) Concentration profiles of sulfate (red) and methane (green) from ODP leg 207, site 1258, Demerara Rise, Equatorial Atlantic 
(Shipboard Scientific party, 2004; Arndt et al. 2006) show a sulfate-methane transition zone in a sediment depth of 150-300 
m. Substrates are supplied by a downward flux of sulfate from the sediment-water interface and by the biogenic production 
of methane from deeply buried organic matter in black shale sequences at a depth below 400 m. Upward diffusing methane 
and downward diffusing sulfate are depleted by deep sedimentary microbial communities in the process of anaerobic methane 
oxidation (AMO). Note that the exponential distribution of methane ranges over several orders of magnitude. 

(b) Depth profiles of sulfate (red) and barium (blue) along the sediment column of ODP leg 201, site 1228, eastern Pacific 
ocean (Shipboard Scientific party, 2004; D'Hondt et al, 2004), exhibiting two reversed zonation patterns. The characteristic 
sulfate-barium transition zone that extends from the water-sediment interface to greater depth (0 to 40 m) is mirrored by a 
second reversed succession that extends upward from the basement-sediment interface in depth from 40 to 80 m. Here, sulfate 
enters the sediment from two directions: from the overlying ocean and from an underlying basaltic aquifer. In each transition 
zone the depletion of sulfate by microbial activity promotes the remobilization of biogenic barium (Torres et al. 1996), giving 
rise to high concentrations of dissolved barium in the pore fluids beyond the zone of sulfate depletion. 

Data are obtained during the ocean drilling program (ODP, http://iodp.tamu.edu/[ ). Dashed and dotted lines show straight 
line fits to the data on the logarithmic axes. values of the fit are indicated. 
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Figure S6: Illustration of different competition outcomes in the phytoplankton model. The outcome of two species competition 
is shown (indicated by color coding), depending on the critical resource values of species 2 (critical resources (A'^i , Ii) of species 
1 fixed, green circle), when the maximal growth rates and mortalities are identical (top panel), or differ by the factor P — 2 
(bottom panel), similar to Fig. 5 from the main text. Insets (a) - (d) show typical spatial profiles corresponding to the different 
competition regimes (critical resource values of species 2 indicated as red circle). Plotted are the phytoplankton concentration 
of species 1 (green) and species 2 (red), and the distribution of light (black dashed line) and nutrient (black solid line) as a 
function of depth z. 

(a) and (d): Two fundamentally different coexistence mechanisms in a spatial system, (a) Coexistence due to a resource 
limitation trade-off, mediated by niche segregation in resource requirements which becomes apparent as a spatial separation 
of density profiles, (d) Coexistence due to a gleaner-opportunist trade-off. The two species are not spatially separated, but 
species 1 with smaller resource requirements (gleaner) can utilize a larger favorable range, whereas the high growth rate of the 
stronger resource limited species 2 (opportunist) allows it to survive on a smaller spatial range. 

(b) and (c) ; Two kinds of bistability in the competition outcome, (b) Alternative stable state when each species cannot grow 
in the presence of its competitor, (c) Alternative stable states of either coexistence (species 1 can invade the monoculture of 
species 2, but does not reach a high abundance) or a monoculture of species 1 (species 2 cannot establish in the presence of 
species 1). 

Parameter values of species 1 see in Table [ST] Simulation parameters for species 2: (a) Hi^2 ~ 24.7 /imol photons m~^ s^^, 
Hn,2 = 0.012 mmol nutrient m~^, (b) Hi^2 = 12.2 ^mol photons m~^ s~^, Hn,2 ~ 0.123 mmol nutrient m~^, (c) Hi^2 ~ 9 /imol 
photons m~^ s~^, Hn,2 ~ 0.22 mmol nutrient m~^, (d) Hi ^2 ~ 24.7 /imol photons m~^ s~^, Hn,2 ~ 0.044 mmol nutrient m~^. 
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Table SI: Default model parameters 



Symbol Interpretation 




Units 


Value 








Independent variables 














t Time 




h 










z Depth 




m 










Dependent variables 


















L-CliB 111 










I{zA) Light intensity 




/imol photons m"^ s^^ 










N{z, t) Nutrient concentration 




mmol nutrient m""^ 










Parameters 














/in Incident light intensity 




/imol photons 


600 (100 


. . . 1000)* 






is'bg Background turbidity 




m-i 


0.1 (0.01 


...1)* 






Zb Depth of the water column 




m 


100 








D Vertical turbulent diffusivity 




cm^ 


0.3 (0.01 


...0.7)* 






Nb Nutrient concentration at Zb 




mmol nutrient m"^ 


10 (0.1.. 


.10)* 






Species 1 














Hi half-saturation constant for lij 


;ht 


/imol photons m"^ s^^ 


20 








Hn half-saturation constant for nutrient 


mmol nutrient m""^ 


0.04 








k hght attenuation coefficient 




m^ ceir^ 


6x10"^° 


(1x10-1° . 


. . 1x10" 




a cell nutrient content 




mmol nutrient celP^ 


2xl0~i° 


(1x10-1° . 


. . 1x10" 




Mmax maximal growth rate 




h-i 


0.04 








m mortality rate 




h-i 


0.01 








Species 2 














Hi half-saturation constant for lif 


;ht 


/imol photons m"^ s"^ 


15 








Hn half-saturation constant for nutrient 


mmol nutrient m~"^ 


0.065 








k light attenuation coefficient 




m^ cell-i 


1x10-1° 


(1x10-1° . 


. . 1x10" 


-9). 


a cell nutrient content 




mmol nutrient cell^^ 


5x10^1° 


(1x10-1° . 


. . 1x10" 




Mmax maximal growth rate 




h-i 


0.04 








m mortality rate 




h-i 


0.01 









* in Fig. 3b and Fig. [Ml 

Appendix S5 provides further details of simulation parameters specific to certain figures. 



